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CHAPTER  I 

INTRODUCTION 

Maintainability  (M) ,  as  an  engineering  discipline,  is 
progressing  through  its  infancy  and  is  now  gathering  sophistication 
and  momentum.  This  growth  can  be  seen  in  the  increasing  number  of 
published  articles  introducing  new  methods,  new  applications  of  old 
methods,  and  new  concepts  slanted  toward  the  field  of  maintainability. 
This  paper  can  be  placed  in  the  second  category.  It  is  an  old 
concept,  applied  in  a  different  way  to  obtain  the  desired  results. 

The  desired  result  in  this  case  is  the  analysis  of  downtime 
data  by  fitting  a  distribution  to  the  data.  The  old  method  is  the 
familiar  gamma  distribution  which  has,  in  recent  years,  been  the 
subject  of  lively  conversations  as  a  possible  distribution  of 
downtimes.  It  has  been  pointed  out  by  many  authors,  Peterson  (14)*, 
in  particular,  that  the  gamma  distribution  can  assume  a  variety  of 
different  shapes.  This  is  verified  by  Figure  1.  In  particular,  note 
tl  at  when  the  shape  parameter  is  equal  to  1,  the  distribution  is  aa 
exponential.  Although  the  shape  parameter  is  not  restricted  to  an 
integer  value,  an  Erlang  distribution  (1)  results  when  it  is.  As 
the  Erlang  distribution  has  boon  used  in  queueing  theory,  ome 
further  attributes  of  the  gamma  distribution  may  be  for  there. 

Peterson  C14) ,  in  his  paper,  discusses  the  fittir  ;  of  a  gamma 
distribution  to  empirical  data  generated  bv  a  computer.  Although 
he  discusses  the  attributes  of  the  gamma  distribution  in  the  analysis 

*  Numbers  in  parentheses  re  lor  to  the  material  listed  in  the 
References  at  the  end  of  this  paper. 
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of  downtime  to  quite  an  extent,  his  method  of  parameter  estimation 
leaves  much  to  be  desired  statistically.  This  method  consisted  of 
plotting  the  empirical  data  and  then  "overlaying"  this  plot  with 
plots  of  the  gamma  distribution  with  integer  shape  parameters.  The 
gamma  distribution  which  fitted  best  by  this  overlaying  process  was 
then  used  to  solve  for  tVj  scale  parameter.  The  threshold  parameter 
was  assumed  to  be  zero  in  all  cases.  Although  the  results  were 
biased  by  this  method,  he  does  give  an  actual  application  of  the 
gamma  distribution  fitted  to  downtime  data. 

The  gamma  distribution  is  attractive  for  two  main  reasons. 

The  first  is  that  the  log-normal  distribution,  which  is  usually 
assumed  to  be  the  distribution  of  downtimes,  is  quite  unwieldy  and 
requires  plotting  on  log-normal  paper  to  ascertain  that  the  data  is 
actually  log-normally  distributed.  The  second  reason  is  that  many  of 
th  tools  and  methods  being  brought  into  use  today  ..ometimes  require 
the  manipulation  of  distributions  in  complex  ways.  For  example, 
Pieruschka  (15)  develops  the  force  of  mortality  of  equipment  with 
ideal  repair  and  shows  that  a  steady  state  value  is  obtained  as  the 
number  of  failures  increase.  He  uses  Laplace  transforms  to  obtain 
the  convolution  of  a  gamma  distribution  with  integer  shape  parameters 
and  then  integrates  the  resulting  distribution  to  find  the  mean. 

These  are  the  reasons  whv  i(  is  desirable  to  prove  that  the  eamma 
distribution  will  represent  the  anticipated  data,  and  still  be 
mathematically  simple. 

Two  downtime'  distributions  can  he  spoken  of  in  equipment  analysis. 
The  first  is  the  distribution  of  elapsed  time  required  to  complete 
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a  specific  maintenance  task..  The  second  concerns  the  distribution 
of  times  required  to  complete  different  tasks  within  a  given  system, 
and  this  distribution  is  said  to  describe  the  system  downtime.  The 
random  variable  in  the  second  c.ise  is,  of  course,  composed  of  elements 
of  both  cases,  but  the  randomness  of  the  individual  task  times  is 
usually  neglected.  This  time  element  approach  to  the  analysis  of 
downtimes  will  not  be  used  here;  but,  instead,  we  will  only  be 
interested  in  fitting  a  function  to  a  number  of  data  points. 

The  function  to  be  fitted  is  the  gamma  distribution  and  para¬ 
meter  estimation  is  by  maximum  liktlihood  where  certain  terms  of  the 
maximum  likelihood  equations  have  been  tabulated.  These  tables  are 
entered  with  a  parameter  formed  by  the  ratio  of  the  arithmetic  mean 
to  the  logrithmic  mean  and  the  shape  parameter  term  is  read  out. 
Several  more  simple  calculations  will  yield  the  threshold  parameter 
and  the  scale  parameter. 

In  obtaining  probability  levels,  the  gamma  distribution  has  one 
of  the  disadvantages  of  the  leg-normal  distribution  in  that  tables 
must  be  consulted  {8,li'.  All  of  the  moments  of  the  gamma  distribu¬ 
tion  are  functions  or  ’-.re-  parameters,  and  possiblv  most 

important,  a  gamir.c.  cisti ibution  ccnvoi  ^d  with  itself,  results  in 
another  gamma  distribution. 

To  prove  that  the  gamma  distribution  cojld  be  fitted  to  downtime 
data,  a  computer  program  was  written  to  perform  the  estimation  pro- 

...  .  •  r-i  : v  -  — :pl  ••  e  iy  specified  density  functions 

*.r,  a.  cm  pi  r  Lc  a  i  *.  r  fit  w  t  h  ;  gamma  distribution.  When 

tested  by  the  Kolmogorov-. -r  rnoff  goodness  of  fit  test  (11),  there 
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was  no  significant  difference  between  any  of  the  sets  of  data  and  its 
corresponding  gamma  distribution  approximation  at  the  .05  level. 

A  brief  literature  survey  will  be  given  in  Chapter  II,  with  the 
ensuing  chapter  devoted  to  the  development  and  use  of  the  gamma 
distribution  and  the  estimation  procedure.  The  convolution  of  the 
gamma  distribution  as  pertains  to  downtime  data  is  developed,  and 
several  examples  will  be  given  in  Chapter  IV.  Some  possible  extensions 
of  and  conclusions  to  the  paper  will  then  be  discussed  in  Chapter  V. 
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CHAPTER  II 

LITERATURE  REVIEW 

Much  of  the  impetus  for  this  paper  resulted  from  Peterson's 
paper  (14)  which  was  a  rough  approximation  of  downtime  data  by  the 
gamma  distribution.  Initially,  work  was  started  on  the  so-called 
generalized  gamma  distribution  developed  by  Stacy  (16)  of  the  IBM 
Corporation  and  at  times  worked  on  by  Stacy  and  Mihram  (17),  Harter 
(7),  Harter  and  Moore  (6),  Parr  and  Webster  (12),  and  Collins  (3). 
The  fourth  parameter  is  a  power  of  the  exponential  term  in  the  three 
parameter  gamma  distribution.  Although  Stacy  and  Mihram  and  Collins 
have  developed  two  different  parameter  estimation  procedures,  the 
methods  are  unwieldy  and  time  consuming.  The  convolution  of  the 
generalized  gamma  distribution  is  given  as  a  summation  term.  In 
attempting  to  fit  the  four  parameter  gamma  distribution  to  data,  it 
was  soon  realized  that  the  advantage  of  a  wider  variety  of  possible 
shapes  was  net  worth  the  add . t i era  1  complexity. 

Although  the  method  of  moments  in  parameter  estimation  for  the 
gamma  di;.t  .on  i-  q,  •  •  i  -  im  b  ,  r k.-  greater  confidence  obtained 

by  maxi. -.urn  like  i  .  4  i  •  ..il  lv  wort  i.  It-.  Several  tables  for 

solving  the  maximum  likelihood  '-quations  are  available.  Among  these 
are  Chapman  (2">,  Masuyamn  and  Kuraiwa  C1951),  and  Greenwood  and 
bn rand  tS),  the  last  being  used  in  this  paper.  General  discussions 
t »  •  :i  .••'iuiti  i-c  contained  in  most  statistics  and  relia- 

,  ‘  ,  :u  ;  la  .i  ,  i  n.i',  d.ivelcpeti  the  gamma  distribution 

*_  •  qu:  •’  .«•;  i  Ater.t  in  tei  .  analysis,  and  develops  the 
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mortality  and  hazard  rates  of  the  distribution.  Probability  tables 
for  the  Incomplete  Gamma  Function  have  been  computed  by  Pearson  (13) 
and  Harter  (8).  Some  of  these  discussions  are  presented  in  the  next 
chapter  in  a  general  discussion  of  the  gamma  distribution. 


CHAPTER  III 


THE  GAMMA  DISTRIBUTION 


In  this  chapter,  the  gamma  distribution  will  be  investigated. 
The  moments  and  the  moment  generating  function  along  with  the 
convolution  of  the  gamma  distribution  will  be  given  and  the  method 
of  parameter  estimation  will  be  developed. 

The  gamma  density  function  is  given  by 


f  (t) 


a  f(b) 


1  0 


c  "  t  <  « 

otherwise 


where  c  is  the  threshold  value, 
a  is  the  scale  pa  tame ter, 
b  is  the  shape  par.muli r,  and 
r(n)  is  the  gamma  function  with  argument  n  (Table  2). 

The  cummulative  distribution  is 

l*(tf,)  •  ir  (t  -  tn)  =  I  (u,p) 

where  T.  ui.pl  is  the  Inc  "'.pletn  Gamma  Function  Ratio  and  in  Harter's 
■notation  1  H  > 

p  =  h  -  1 

a.'  '!>-! 

■.*  it  •  ■  > .  :...ik  ■  a !  :  i’  >unJ  t  e  e. 


mortality  and  hazard  rates  of  the  distribution.  Probability  tables 
for  the  Incomplete  Gamna  Function  have  been  computed  by  Pearson  (13) 
and  Harter  (8).  Some  of  these  discussions  are  preset  ed  in  the  next 
chapter  in  a  general  discussion  of  the  gamma  distribution. 
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•  (b  +  k). 
'  (b) 


k-1 

ak  ::  (h  +  j). 

i=o 


A  point  that  will  bo  used  later  should  be  noted  here,  and  that  is 
that  a  number  of  different  combinations  of  scale  parameter  and  shape 
parameter  will  yield  approximately  the  same  distribution.  This  can 
be  seen  from  the  equations  for  the  mean  and  variance  in  which  there 
are  three  unknowns.  An  illustration  of  this  point  is  shown  in 
Figure  2  where  the  respective  parameters  of  the  curves  are 


Curve 

c 

a 

b 

Mean 

Variance 

1 

- .  2 

.65 

A.  15 

2.5 

1.75 

2 

0 

.707 

3.53 

2.5 

1.76 

3 

# 

.89 

2.36 

2.5 

1.87 

Parameter  Kstimation  hv  Maximum  Likelihood 


In  the  development  ot  the  maximum  likelihood  equations  for 
parameter  estimation  fV),  we  form  the  joint  probability  of  drawing 
n  samples,  tj,  t7,  t-^...,  tn  I  ram  a  gamma  distribution  with  parameters 
c,  a,  and  b. 


FIGURE 
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n 

n 

1=1 


(ti-c) 


b-1 


e 


n 
-  I 
i  =  l 


ti~c 

a 


ab"  r_r(b)ln 


Taking  the  logarithm  yields 

n  n  t 

InL  =  (b-1)  ln(t  -c)  -  T.  — — —  -  bn  In  a  -  n  F  (b) . 

i=l  i=l  a 

Differentiating  with  respect  to  the  parameters  to  maximize  the 
probability  of  obtaining  the  given  sample  results  in 

=  ti~c  _  nb  _  0 
*a  1=1  a  '  d 

,  I  t-f  -  c 

ab  =  - - -  Equation  1 

=  Z  In  ft.  -  c)  -  n  In  a  -  n  i — =  0. 

^b  i=l  :'h) 

Substituting  the  value  found  for  a  in  equation  1,  one  obtains 

n 

.  In  (t^-c) 

In  b  -  1  1 (b)  =  In  ^i~C)  -  Izi _  Equation  2 

F  (b)  n  n 

where  I  ' (b)  is  the  first  derivative  of  the  gamma  function  with  respect 
to  b.  Taking  the  partial  with  respect  to  the  threshold  parameter  gives 

_  =  -  (b  -  1)  v  - 1 —  +  _  =  0 

>C  tf  -  c 


a  (b-1)  - - —  =  n.  Equation  3 

i=l  ci  '  c 


In  equation  1  the  product  ab  is  seen  to  be  the  arithmetic  mean 


corrected  for  the  threshold  value.  In  equation  2,£r'(b)J  *  £7  (b)J 
Is  the  tabulated  digarama  function  (4),  and  equation  2  may  be  solved 
using  a  °et  of  these  tables  and  natural  logarithms.  However,  the 
entire  left-hand  side  of  equation  2  has  been  tabulated  in  several 
different  ways  (2,5). 

This  then  will  be  the  method  used  to  estimate  the  parameters. 

The  right-hand  side  of  equation  2  is  simply  the  logarithm  of  the 
arithmetic  mean  minus  the  logarithm  of  the  geometric  mean.  If  we 
denote  this  right-hand  side  r s  Y  and  the  left-hand  side  as  n(b) ,  then 

n(b)  =  Y. 

Greenwood  and  Durand  (5)  have  tabulated  brt(b)  =  b[ln  b  -  r'(b)/ 

F (b) 7  versus  n(b)  =  Y  to  eight  decimal  places  and  their  table  is 
reproduced  as  Table  1. 

The  procedure  will  be  to  calculate  the  arithmetic  and  geometric 
means  of  the  data  and  set  Y  equal  to  the  difference  of  the  logarithms 
of  these  two  means.  Enter  Table  1  with  this  value  and  read  bn(b)  =  Yb. 
Then  solve  for  b  =  Yb/Y. 

The  scale  parameter  a  is  found  from  equation  1  and  is  equal  to 
the  arithmetic  mean  divid’d  by  b.  Equation  3  indicates  that  some 
iteration  process  must  be  used  to  solve  for  the  threshold  value. 

Since  this  value  must  be  found  before  the  means  can  be  calculated,  a 
method  of  estimation  which  is  free  of  terms  containing  a  and  b  must 
be  found.  The  following  estimate  was  developed  by  fitting  a  straight 
•'  -ie  to  the  firs*-  three  p-.irts  on  the  curve  generated  by  the  expected 
ranks  of  the  various  ordered  sample  points.  (Appendix  3  shows  the 
development  of  this  equation').  The  resulting  expression  is 
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c  -  Y2  +  Y1  .  (Y3  -  Yl>  <Y3  -  Y2> 

2  2(2Y2  -  Y !  -  Y3) 

where  Y^,  Y2,  and  Y2  are  the  first  three  ordered  sample  points. 

Perhaps  the  most  valuable  asset  of  the  gamma  distribution  is  the 
ability  to  convolve  one  gamma  distribution  with  another.  The  easiest 
proof  of  this  is  to  derive  the  moment  generating  function  and  then 
find  the  distribution  of  the  sum  of  random  variables  by  multiplying 
their  respective  moment  generating  functions  together.  For  the  gamma 
distribution,  the  moment  generating  function  is 


Thus,  a  gamma  distribution  convolved  witi  itself  will  result  in 
a  gamma  distribution  with  threshold  parameter  2c,  scale  parameter  a,  and 
shape  parameter  2b.  In  general,  when  n  gamma  distributions  are  con¬ 
volved  together,  if  they  have  a  common  scale  parameter,  a  gamma  distri- 

n 

butic.i  results  ^ith  threshold  parameter  =  c.,  scale  parameter  a, 

n  i  =  l 

and  shape  parameter  =  I  b^.  The  restriction  of  equal  scale  para- 

i=  1 

meters  for  the  distributions  being  convolved  will  not  be  too  great  of 
an  impairment  in  the  treatment  of  downtimes  since  each  set  of  data 
will  have  approximately  the  same  range.  Recalling  that  a  set  of  data 
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can  be  represented  by  a  variety  of  scale  and  shape  parameter  combina¬ 
tions,  it  fs  seen  that  sets  of  data  with  unequal  scale  parameters  can 
be  convolved  by  preassigning  a  common  scale  parameter  to  each.  The 
shape  parameter  for  each  set  of  data  would  then  be  found  by  dividing 
thp  individual  arithmetic  means  by  the  common  scale  parameter.  This 
concept  will  be  discussed  in  the  following  chapter  in  the  form  of  an 
example. 
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CHAPTER  IV 

EXAMPLES 


The  first  example  in  this  chapter  shows  the  steps  necessary  to 
fit  a  gairana  distribution  to  a  set  of  data.  The  second  example  will 
illustrate  the  method  of  convolving  several  sets  of  data  with  unequal 
scale  parameters. 

1.  The  lollowing  times  were  taken  in  a  maintenance  test  on  a 
particular  system.  We  wish  to  perform  an  analysis  on  the  data  and 
fit  a  distribution  to  it  so  that  we  may  speak  of  the  probability  of 
completing  a  future  action  or  actions  within  a  certain  time. 


1.7 

4.6 

2.3 

3.0 

4.0 

3.4 

1.3 

3.7 

3.  3 

ON 

• 

CM 

3.8 

2.2 

2.6 

o 

• 

CN 

2.7 

4.3 

3.5 

2.5 

3.1 

The  first  step  in 

the  analysis 

is  to  find  the 

first  three  ordered 

points.  They  are 

Y!  =  1.3,  Y2 

=  1.7,  Y3  =  2. 

0,  then 

c  - 

1.7  +  1.? 

(2.0  -  1.3)  (2 

■0  -  1-7) 

2  2[(2)  (1.7)  -  1.3  -  2.01 


=  1.5  -  -til  =  1.5  -  1.05  =  .45 

.  2 

The  arithmetic  mean  is  found  to  be 


48.35 


x 


n 
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2.542. 


where?  .301  Is  the  critical  val.io. 


GAMMA  DISTRIBUTION  FITTED  TO  NORMAL  DATA 
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Therefore  there  Is  no  significant  difference  at  the  .05  level  of 
significance. 

2.  A  second  example  will  demonstrate  the  convolution  of  a  number  of 
sets  of  data.  The  data  is  as  follows: 

System  1: 


3.13 

1.96 

3.84 

2.33 

2.48 

3.00 

1.36 

3.38 

4.64 

4.04 

2.62 

1.72 

2.87 

2.16 

3.67 

4.28 

3.52 

3.25 

2.75 

System  2: 


4.57 

4.04 

5.76 

6.87 

5.59 

3.81 

4.71 

5.43 

3.54 

5.00 

5.29 

4.24 

6.19 

6.46 

3.13 

5.96 

4.85 

5.15 

4.41 

System  3: 


3.11 

1.42 

1.67 

2.22 

2.45 

2.11 

.58 

.89 

1.89 

1.78 

2.58 

2.90 

3.42 

2.73 

2.00 

1.55 

1.27 

2.33 

1.10 

7.61 

5.16 

4.39 

5.52 

5.35 

© 

o 

7.50 

6.48 

4.50 

7.06 

6.16 

5.68 

6.65 

5.84 

8.07 

6.32 

i.ii 

o  84 

Svstem  4; 


Fitting  the  gamma  distribution  to  the  four  systems  as  in  the  first 


example , 

the  various  parameters  are 

found  to  be. 

System 

c  a 

b 

Mean 

Variance 

1 

.94  .46 

4.5 

3.0 

.95 

2 

2.68  .14 

4.3 

5.0 

1.26 

3 

2.03  .34 

11.7 

6.0 

1.35 

4 

.19  .39 

4.7 

2.0 

.715 

4 

Setting 

i~l  A- 

A  =  1  =  .43, 

4 

the  bj's 

are  then 

bl 

o 

00 

it 

«-*)< 

'*1 

ii 

b2 

=  ^  =■  1.41, 

A 

b3 

=  —  =  9.24,  and 

A 

b4 

=  =  4.21. 

A 

Then  the 

distribution  of  Liie  convolveu  systems  i; 

s  a  gamma 

distribution 

4 

with  threshold  parameter  C  =  c.  =  1.84,  scale  parameter  A  =  .43, 


and  shape  parameter  B  =  =  23.7.  This  distribution  has  been 

i=l 

plotted  in  Figure  4  along  with  the  convolution  ol  the  four  data 
generating  distributions  which  were  all  normal  with  means  of  3,  5,  6  and 
2  and  variances  of  1,  1.3,  1.6  and  7.1. 


Note  in  particular  the  parameters  of  the  total  distribution 
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Method 

Mean 

Variance 

t95% 

fc99% 

Convolved 

C  amma 

16 

4.4 

19.74 

21.59 

Convolved 

Normal 

16 

4.65 

19.56 

21.03 

The  Kolmogorov-Smirnof f  goodness  of  fit  was  again  applied  with 
4  x  19  =  76  degrees  of  freedom.  No  significant  difference  could  be 
found  at  the  ,0j  level. 
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CHAPTER  V 

CONCLUSIONS 

Although  the  gamma  distribution  approximations  found  by  the 
author  fitted  extremely  well,  it  is  recommended  that  someone  using 
the  method  go  slowly  and  check  the  initial  results  before  proceeding. 
This  method  will  have  important  ramifications  in  the  design  of  new 
systems  since  the  designer  will  know  beforehand  that  the  downtimes  will 
follow  a  gamma  distribution  of  some  sort.  He  can,  therefore,  feel  free 
to  use  the  models  necessary  for  proper  design  of  a  logistical  system 
without  making  a  poorly  based  assumption  as  to  the  distribution  of  the 
downtimes.  Also,  with  the  many  time  critical  repair  action  situations 
that  occur  in  the  space  industry,  a  planning  engineer  can  use  the 
accurate  probability  models  without  employing  an  inaccurate  repair 
time  distribution. 

The  convolution  tecnniqiu  developed  for  unequal  scale  parameters 
needs  more  validation  than  the  single  example  presented.  However,  the 
potential  of  such  a  technique  demands  that,  it  he  investigated  further. 

For  those  who  would  still  argue  the  validity  of  the  log-normal 
distribution  in  maintainability  work.  Figure  U  has  been  provided.  It 
is  a  gamma  distribution  fitted  to  a  log-normal  distribution  with  a 
moan  of  ,7h  and  variance  of  .604. 

Tn  the  estimation  procedure,  a  and  b  were  fixed  hv  the  data  as 
soon  as  the  threshold  parameter  was  found,  and  varying  c  will  vary 
:he  fit.  The  equation  for  c  on  page  11  is  certainly  not  the  best 
estimate.  A  possible  alternative  would  be  to  plot  the  first  few 


Gamma  (.5,  .966,  2.07) 


Distribution  Fitted  to  Log-Nonnal 


sample  points  versus  their  expected,  or  median,  ranks  and  approximating 
the  threshold  parameter  by  eye.  This  would  certainly  lend  self-assurance 
to  the  estimate. 


Although  the  process  can  be  hand  operated,  computer  programs  have 
been  written  for  both  the  data  fitting  process,  and  for  the  convolution 
technique  and  are  presented  with  the  expected  output  in  Appendices  1 
and  2. 

Greenwood  and  Durand  C5)  have  also  developed  the  formulas  for  the 
variances  of  a  and  b.  Since  we  violated  one  of  their  original  premises 
by  estimating  c  from  the  data,  these  formulas  will  not  be  presented 
in  order  to  avoid  temptation  of  the  reader. 

Further  areas  of  research  in  the  line  of  using  the  gamma  distribu¬ 
tion  for  representing  repair  times  would  be  to  improve  the  estimate  of  the 
threshold  parameter  in  order  to  obtain  a  better  fit.  Perhaps  this  should 
be  developed  for  a  better  fit  in  the  upper  half  or  perhaps  only  the  upper 
tail.  The  convolution  technique  has  assumed  equal  weight  for  all  of  the 
systems,  this  could  most  certainly  not  be  the  case.  Further  research 
could  be  done  on  weighting  the  different  sets  of  data  by  their  frequency 
with  respect  to  the  overall  system.  Interesting  possibilities  also  lie 
in  the  previous  research  done  on  the  Erlang  distribution  (1)  which  is  a 
special  case  of  the  gamma  distribution.  Perhaps  some  of  the  theories 
which  have  been  developed  in  queueing  theory  about  the  Erlang  distribu¬ 
tion  would  be  useful  in  maintenance  queueing  analysis. 

It  is  hoped  bv  the  author  that  further  research  on  the  gamma  distri¬ 
bution  will  be  done  and  that  the  gamma  distribution  will  be  used  more 
frequently  in  downtime  analvsis.  Since  the  gamma  distribution  is  quite 


Ji 
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simple  in  form,  especially  when  the  shape  parameter  is  an  integer,  it 
is  easy  to  become  familiar  with  the  distribution  and  gain  an  insight 
into  the  data  it  represents.  It  is  hoped  that  this  paper  will  help 
make  the  gamma  distribution  an  available  tool  of  the  Maintainability 
engineer. 


Y  VERSUS  Yb 
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APPENDIX  1 

Computer  Program  for  Fitting  a  Gamma 
Distribution  to  a  Set  of  Data 

This  appendix  contains  the  flowchart,  the  program,  and  the 
expected  output  of  a  program  which  fits  a  gamma  distribution  to  n 
data  points.  In  addition,  the  input  formats  are  given  on  standard 
forms.  The  process  is  as  in  the  main  portion  of  this  paper  with  the 
following  exceptions. 

1.  In  lieu  of  Tabic  1,  the  following  equations  (5)  were  used  to 
solve  for  b. 

b  =  y_1(. 5000876  +  . I648852y  -  .0544274y2) 

For  0  <  y  .  5772 

b  =  8*898919  +  9.059950V  +  ,9775373v2 
y( 17. 79728  +  11.968477y  +  y^) 

For  .5772  --  y  <  17.0 

2.  The  geometric  mean  was  calculated  by 


FLOWCHART 


o  n  n  n 


DIMENSION  X(100),Z(100) 

RE AD (5, 10)  N 

READ (5, 20)  (X (M) ,  M=1,N) 

C  THERE  ARE  N  SAMPLE  JOINTS  X(l) ,X(2) , . . .X(N) 

SXBAR  =0.0 
SXDOT  =1.0 
AN  =  N 

C=(X(2)+X(l))/2.0-(X(3)-X(l))*(X(3)-X(2))/(4.*X(2)-2. 
1  *X(1)-2.*X(3)) 

C  C  IS  THE  THRESHOLD  PARAMETER 

DO  100  L  =  1,  N 
Z (L)  =  X(L) -C 
SXBAR  -  SXBAR  +  Z(L) 

SXDOT  =  SXDOT  *  Z(L)**(1.0/AN) 

100  CONTINUE 

XBAR  =  SXBAR/AN 
Y  =  ALOG (XBAR/ SXDOT) 

C  XBAR  IS  THE  ARITHMETIC  MEAN,  SXDOT  THE  GEOMETRIC  MEAN 

IF (Y-. 5772)  1,1,2 

1  B  =  (.5000876+. 1648852*Y- ,0544274*Y*Y)/Y 
GO  TO  3 

2  B  =  (8. 898919+9. 05995*Y+.9775373*Y*Y) / (Y*<17. 79728 
1  +  11.968477  *  Y  +  Y*Y)) 

C  B  IS  THE  SHAPE  PARAMETER 

3  A  =  XBAR/B 
A  IS  THE  SCALE  PARAMETER 

VALUE  OF  THE  DENSITY  FUNCTION  AT  EACH  SAMPLE  POINT. 
THE  NEXT  SIXTEEN  CARDS  ARE  OPTIONAL  AND  WRITE  THE 
THE  FIRST  SEVEN  CARDS  FIND  GAMMA 'B) 

WRITE (6, 30) 

30  FORMAT (1H1,14X, 'SAMPLE ' ,GX, 'SAMPLE ’ , 13X , ’DENSITY' / 

1  14X, ’NUMBER’ ,9X, 'POINT'//) 

AMUL  =  1.0 
BG=B 

7  IFfBG  -  2.0)  4,4,5 

5  BG  =  BG  -  1.0 

AMUL  =  AMUL  *  BG 
GO  TO  7 

4  Q  =  BG  -  1.0 

GAMB=AMUL*  G 1 .  - .  57710 166*Q+.  98 58 5  399*Q**2 .  - .  87642 128*Q 

1  **3 .+ . 8  328  2 1 2*Q**4 . - . 5684  729*Q** 5 .+ . 2  548  204  9*Q**6 . - 

2  .0514993*Q**7.) 

DO  200  K  =  1,N 

FX  =  ((ZGK)/A)**(B-1.)*EXPGG-1.)*Z(K)/A))/(A*GAMB) 
WRITE ( 6 ,40)  K,  X(K),  EX 
40  FORMAT  <\ 5X , 1 3 , 10X , E6 . 2 , 5X , F 16 . 7 ) 

200  CONTINUE 

WRITE  <  6 , 50)  C, ,  A,  B 
10  FORMAT  HI) 

20  FORMAT (12F6. 2) 


34 


50  FORMAT (1H0, 'THE  DATA  IS  GAMMA  DISTRIBUTED  WITH  ' , 

1 'THRESHOLD  PARAMETER  =  ' .F7.4//40X, 'SCALE  PARAMETER  -  ' 
2  .F8.4//36X, 'AND  SHAPE  PARAMETER  =  ',F8.4) 

STOP 

END 
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CHAPTER  III 

THE  GAMMA  DISTRIBUTION 


In  this  chapter,  the  gamma  distribution  will  be  investigated. 
The  moments  and  the  moment  generating  function  along  with  the 
convolution  of  the  gamma  distribution  will  be  given  and  the  method 
of  parameter  estimation  will  be  developed. 

The  gamma  density  function  is  given  by 


f  (t) 


a  f(b) 


0 


c  <  t  <  *° 


otherwise 


where  c  is  the  threshold  value, 
a  is  the  scale  parameter, 
b  is  the  shape  parameter,  and 
;(n)  is  the  gamma  function  with  argument  n  (Table  2). 

The  eiiinmululive  distribution  is 

F(tr)  «  Pr  (t  t0)  =  I  (u.p) 

where  T.  <  u  M  is  the  Incomplete  Gamma  Function  Ratio  and  in  Harter's 
notation  <  “■ 1 

p  "  h  -  1 

•  i  -  ^  °  1 

a*  h-1 

a t  i  >unJ  :  ,  iv 


■  i: 


* »  :  t 


0.0659849 


NPUT  FORMAT  COLUMN  IDENTIFICATION 


NUMBER  OF  SAMPLE.  POINTS,  N,  FORMAT  (13)  0  19 


NUMBER  OF  SAMITE  POINTS 


38 

APPENDIX  2 

Computer  Program  for  Convolving 
Several  Sets  of  Data 

This  program  utilizes  the  previous  program  in  a  do  loop  to  fit 
gamma  distributions  to  several  sets  of  data.  The  average  scale  para¬ 
meter  is  found  and  each  set  of  data  is  forced  to  have  this  common 
scale  parameter.  The  distributions  are  convolved  with  the  threshold 
k  k 

parameter  being  E  ci,  scale  parameter  =  I  and  shape  parameter  = 

i*l  i=l 

k 

k 

Z  where  there  are  k  sets  of  data. 
i=l 


Read  the  Number  of  Subsystems,  K 
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DIMENSION  X(IOO) 

READ (5, 10)  K 

C  THERE  ARE  K  SUB-SYSTEMS  BEING  ANALYZED 

SUMC  =  0.0 
SUMA  =  0.0 
SUMB  =  0.0 

C  THE  DO  800  LOOP  FITS  A  GAMMA  DISTRIBUTION  TO  EACH  SET 

C  OF  DATA. 

DO  800  I  =  1,K 

READ (5, 10)  N 

READ ( 5 , 20)  (X (M) ,  M=1,N) 

C  THERE  ARE  N  SAMPLE  POINTS  X(l) ,X(2) , . . „X(N) 

SXBAR  =0.0 
SXDOT  -  1.0 
AN  =  N 

C=(X(2)+X(l))/2.0-(X(3)-X(l))*)X(3)-X(2))/(4.*X(2) 

1  -  2.*X(1)  -  2,*X(3) ) 

C  C  IS  THE  THRESHOLD  PARAMETER 

DO  100  L  =  1,  N 

Z(L)  =  X(L)-C 
SXBAR  =  SXBAR  +  Z (L) 

SXDOT  =  SXDOT  *  Z(l.)**(  1.0/AN) 

100  "ONTINUE 

XBAF  =  SXBAR /AN 
XB(I)  =  XBAR 
Y  =  ALOO( XBAR /SXDOT) 

C  XBAR  IS  THE  ARITHMETIC  MEAN,  SXDOT  THE  GEOMETRIC  MEAN 

IFFY-. 5772)  1,1,2 

1  B  =  ( .  5000876+ .  16488  5  2-'Y-  . 0544274*Y*Y)  /Y 
GO  TO  3 

2  B  =  (8. 898919+9. 05993--Y+.  9775373*Y*Y)  /  (Y*  ( 17 . 79728 
1  +  11. 96 8 A 7 7  *  Y  +  Y*Y) ) 

C  B  IS  TIB:  SHAPE  PARAMETER 

3  A  =  XBAK/B 

C  A  IS  THE  SCALE  PARAMETER 

SUMA  =  SirMA  +  A 
SUMC  =  SUMC  +  C 

C  THIS  WRITE  STATEMENT  IS  OPTIONAL  AND  PRINTS  THE 

C  PARAMETERS  OF  EACH  SUB-SYSTEM. 

WRITE (6, 30)  1,  G,  A,  B 

30  FORMATdHO, ’SUB-SYSTEM' ,13, MIAS  BEEN  FITTED  WITH  A', 

1'  GAMMA  DISTRIBUTION  AND'//'  G  ='F8.4,'  ,  A  ='F8.4, 

2'  AND  B  =  ' .F8.41 
800  CONTINUE 

AK  =  K 

AAVG  =  SUMA/AK 
DO  900  L  =  1 , K 
SUMB  =  SUMB  +  XB(L)/AAVG 
900  CONTINUE 

WRITE (6,40)  K,  SUMC,  AAVC,  SUMB 
10  FORMAT (1 3) 
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20 

40 


FORMAT (12F6. 2) 

FORMAT (1H°, I 3, '  SUBSYSTEMS  HAVE  BEEN  FITTED 
nJ^^BUTI0NS  AND  THEIR  CONVOLUTION  IS’// 
3F8  ^STRIBUTI0N  WITH  c  -  ' ,F8.4, '  A  =’,F8.43 
STOP 


WITH  GAMMA’ 
'  A  GAMMA' 

’  AND  B  =' 


SUB-SYSTEM  1  HAS  BEEN  FITTED  WITH  A  GAMMA  DISTRIBUTION  AND 
C  =  0.9400  .  A  =  0.4600  AND  B  =  4.4778 


IN’ PUT  FORMAT  COLUMN  IDENTIFICATION 
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•1  BE R  OF  SUBSYSTEMS 


INPUT  FORMAT  COLUMN  IDENTIFICATION 


NUMBER  OF  SAMPIJ.  POINTS 


SAMPIP.  VALUES,  X  ( I ) 
FORMAT  (12F6.2)  1.36 
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APPENDIX  3 

Derivation  of  the  Equation 
for  the  Threshold  Parameter 

As  indicated  in  Chapter  III,  the  maximum  likelihood  equation  for 
the  threshold  parameter  requires  an  iteration  process.  As  this  para¬ 
meter  is  rarelv  known  in  maintainability  analysis,  some  estimate  of 
this  parameter  must  be  made  from  the  data  available.  Since  the  area 
between  adjacent  ordered  sample  points  is  known  fr^m  order  statistics 

to  be  1  CIO),  a  straight  line  can  be  fitted  to  the  unknown  curve 
n+1 

and  the  threshold  parameter  estimated. 
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The  area  between  the  first  and  second  sample  points  is 

— ~  -  s^2  -  yi) 

n  +■  i 

and  the  area  between  the  second  and  third  points  is 


— L_  =  r(y3  -  y2). 
n  +  1 


The  slope  of  the  line  is  then 


slope  =  8 ) 

y3  -  yi 


and  at  the  point  S 


s  =  A  +  slope 


y2  +  yi 


a  ,  y2  +  yi 

A  =  s  -  slope  —4— — i 


Equation  4 


At  the  point  c 

0  =  A  +  slope  (c) 
A  =  -  slope  (c). 

Equating  equation  4  to  equation  5 


Equation  5 


s  -  slope 


y2  +  yl  _ 


=  -  slope  (c) 


_  >’2  +  yi  _ 


slope 


>'2  f  yi  _  (y3  -  y2>  (>’?  -  yp 

l  2(2y2  -  -  y3) 


and 
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CHAPTER  III 

THE  GAMMA  DISTRIBUTION 

In  this  chapter,  the  gamma  distribution  will  be  investigated. 

The  moments  and  the  moment  generating  function  along  with  the 

convolution  of  the  gamma  distribution  will  be  given  and  the  method 

of  parameter  estimation  will  be  developed. 

The  gamma  density  function  is  given  by 

D-l  _  t~C 
e  a 

- - - c  <  t  <  *° 

a  T(b) 

1  0  otherwise 


where  c  is  the  threshold  value, 
a  is  the  scale  parameter, 
b  is  the  shape  paraniL-ler.  and 
p(n)  is  the  gamma  function  with  argument  n  (Table  2). 

The  cummulative  distribution  is 

F(tc)  =  Tr  (t  <  tG)  =  I  (u,p) 

where  I  is  the  Incomplete  Gamma  Function  Ratio  and  in  Harter's 

notation  1 S i 

p  =  b  -  1 

_  |  _  t  o  ~  . 

n»  b-1 

>:  1  me  -  a:  •  found  !  o  i'e 


.*  lie 


